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Dirac electrons in finite graphene samples with zigzag edges under high magnetic fields (in the 
regime of Landau-level formation) are investigated with regard to their bulk-type and edge-type 
character. We employ tight-binding calculations on finite graphene flakes (with various shapes) to 
determine the sublattice components of the electron density in conjunction with analytic expressions 
(via the parabolic cylinder functions) of the relativistic-electron spinors that solve the continuous 
Dirac- Weyl equation for a semi-infinite graphene plane. Away from the sample edge, the higher Lan- 
dau levels are found to comprise exclusively electrons of bulk- type character (for both sublattices); 
near the sample edge, these electrons are described by edge- type states similar to those familiar from 
the theory of the integer quantum Hall effect for nonrelativistic electrons. In contrast, the lowest 
(zero) Landau level contains relativistic Dirac electrons of a mixed bulk-edge character without an 
analog in the nonrelativistic case. It is shown that such mixed bulk-edge states maintain also in 
the case of a square flake with combined zigzag and armchair edges. Implications for the many- 
body correlated-electron behavior (relating to the fractional quantum Hall effect) in flnite graphene 
samples are discussed. 

PACS numbers: 71.70.Di,73.22.Pr, 73. 21. La, 73.43.Cd 



I. INTRODUCTION 

In the last few years, following the isolation of 
monolayer^ and the fabrication of epitaxial^ graphene, 
the physical properties of graphene nanostructures (in- 
cluding elongated graphene ribbons and finite graphene 
samples and flakes) have established themselves as a 
major research direction in condensed-matter physics. 
This development was propelled by theoretical predic- 
tions (see, e.g. Refs. y-Q) that the electronic properties 
of graphene nanostructures are strongly affected by the 
presence and termination character (in particular, zigzag 
or armchair) of the graphene edges, suggesting an unpar- 
alleled versatility and potential for future nanoelectronics 
applications. Crucial to the realization of this perspec- 
tive is the capability to characterize and engineer edges 
with high purity and perfection, a need that has spurred 
an ever expanding experimental effort which has already 
yielded highly promising resultSjiSrJ^ 

In this context, a recent study^^ of ours addressed the 
influence of graphene edges upon the properties of corre- 
lated many-body fractional-quantum- Hall-effect (FQHE) 
states of Dirac electrons. The most recent experimental 
observatioi*i^"— of such FQHE correlated states [in sus- 
pended monolaye r ^^d^ and bilayer^'5 graphene samples 
under high magnetic fields (B)] has marked another mile- 
stone in demonstrating the potential of graphene not only 
for future technological applications, but also for study- 
ing novel fundamental physics behavior. In particular, in 
Ref. Xli we showed that the Dirac-electron spinors in the 
lowest Landau level display a mixed bulk-edge character 
for graphene samples with zigzag edges, with the bulk 
component giving rise to the ly — 1/3 FQHE state (but 
with an attenuated strength), while the edge component 
is responsible for the insulating behavior observe d^^d^i^^ 



at the Dirac neutrality point. 

Naturally, the complexity of the computational many- 
body treatment of the interelectron repulsion necessi- 
tated the use in Ref. flTi of certain simplified assumptions, 
i.e., a circular shape for the graphene sample and an un- 
interrupted zigzag edge. In this paper, motivated by the 
widespread and ongoing experimental activity on perfect- 
graphene-edge engineering (see above), we present a sys- 
tematic study of the properties of Dirac-electron states 
(with respect to both their bulk-type and edge-type char- 
acter) forming the Landau levels in graphene nanos- 
tructures with more realistic shapes^^r^^l (namely, flakes 
with triangular, hexagonal, and square shapes). To this 
end, we utilize a combination of tight-binding calcula- 
tions on graphene flakes with analytic expressions (via 
the parabolic cylinder function a^^'^^ ) for the relativistic- 
electron spinors associated with the continuous Dirac- 
Weyl equation of a semi-infinite graphene plane. 

We demonstrate that, away from the graphene-flake 
edge, the higher Landau levels contain exclusively elec- 
trons of a bulk- type character (for both sublattices) ; near 
the graphene-fiake edge, these electrons are described by 
edge-type states reminiscent of those familiar from the 
theory of the integer quantum Hall effect for nonrelativis- 
tic electrons <^ In contrast, the lowest (zero) Landau level 
of the graphene flakes contains relativistic Dirac electrons 
of a mixed bulk-edge character without an analog in the 
nonrelativistic case. It is shown that such mixed bulk- 
edge states maintain also in the case of a square fiake 
with combined zigzag and armchair edges. 

The paper is organized as follows: 

Sec. iniis devoted to the description of the methodolo- 
gies employed. Specifically, Sec. Ill Al derives the analytic 
expressions for the Dirac- Weyl spinors in the case of a 
semi-infinite graphene plane with zigzag edge termina- 
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FIG. 1. (Color online) The eigenenergies (specifically the in- 
dex v) that solve the transcendental equation as a function of 
Xc (a) for the K valley [see Eq. and (b) for the K' valley 
[see Eq. (|15|) ]. The almost flat segments of the curves corre- 
spond to Landau levels with energies e « y/2n, n = 0, 1, 2, . . .; 
the approximate symbol « signifies that v does not take inte- 
ger values, but comes extremely close to them. The rising-in- 
energy branches correspond to double-edge states (i.e., states 
of edge character on both the A and B sublattices) ; they cross 
the vertical axis at a;c = for (a) odd integer values and (b) 
even integer values. Note that the K valley (a) exhibits a dis- 
persive (varying with Xc = QvIb) quasiflat band with z/ ~ 0, 
while the K' valley (b) exhibits a dispersionless flat band with 
u = [see text and thick line (magenta color online)]; this is 
the only case when v takes an integer {v = 0) value. 



tion; the solutions for both the K (Sec. Ill A II) and K' 
(Sec. Ill A2p graphene valleys are given. An outline of the 
tight-binding approach used here is given in Sec. Ill Bl 

Our tight-binding results at high magnetic field (con- 
cerning the electron-density components of the two 
graphene sublattices and their interpretation through 
comparison with the continuum-model Dirac-Weyl 
spinors) are presented in Sec. IIII Al for triangular flakes. 
Sec. IIII Bl for hexagonal flakes, and Sec. IIII CI for square 
flakes. 

Finally, Sec. |IV] offers a Summary. 



II. METHODOLOGY 

A. Solutions of the Dirac-Weyl equation for a 
semi-infinite graphene plane 

1. K valley 

For a semi-infinite graphene plane under a perpendic- 
ular magnetic field B (with the graphene plane extend- 
ing for < a: < oo and exhibiting a zigzag edge along 
the y axis at a; = 0), the Dirac-electron wave function 
corresponding to the K valley can be written as a two- 
component spinor (the zigzag boundary condition does 



not couple the two graphene valleys) 

V2 \Xb{x) ) 

In Eq. Qy = ky — Ky, with ky, Ky being the linear 
momenta of the electron and the K valley along the y 
direction; the magnetic length Ib — \/hc/eB. 

With the introduction of reduced (dimensionless) vari- 
ables x/Ib X and Xc — QylB, the continuous Dirac- 
Weyl equation coupling the XAix) and xb{x) compo- 
nents is given by 

-^XB + ix~ Xc)xB = exA (la) 

-^XA- [x- Xc)XA = -£XB, (lb) 
ax 

where the reduced energy e = E/{hvF/lB), with vp 
being the Fermi velocity of graphene. 

In general, and prior to invoking any boundary con- 
ditions (that is considering the complete graphene sheet 
for — cx) < x < cx)), the solutions of the system of coupled 
equations in Eq. ([T]) fall into two classes, i.e., for e 7^ 
and e = 0. 

Solutions for e 7^ 0. In this case, one can multiply 
both sides of Eq. (fTb| with e, and then use Eq. ((Taj) to 
eliminate xa- The result is the following second-order 
equation for xb- 

^xBio + l^^+l-^e^xBio^o, (2) 

where 

C = V2{x - Xc) and ly = (3) 

Eq. ([2]) has the standard form of a Weber differen- 
tial equation, and thus its solutions coincide with the 
parabolic cylinder functions ^^iSSii^ i.e., 

XBiO^CMO, (4) 

where Ci, is a normalization constant. 
Using Eq. (jlal) . the recurrence relation 

and Eq. ([3]), the corresponding A component is given by: 

xa{0 = av^i?.-i(0. (6) 

When i/ is a nonnegative integer, n > 0, the parabolic 
cylinder functions reduce to the familiar wave functions 
of the one-dimensional unconfined harmonic oscillator, 

D^iO = 2-"/2e-«^/^i/„ , (7) 

where are Hermite polynomials. The wave functions 
in Eq. ([7]) have the property D„(±C!o) = 0, appropriate 
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FIG. 2. A bulk-bulk state in the third (n = 2) LL: The xb 
(top frame) and xa (bottom frame) Dirac-spinor components 
are displayed for Xc = 7 and u = 2 + 9.229 x 10"^* [see 
the transcendental Eq. ((8|]. This state lies well inside the 
quasiflat segment of the third LL curve in Fig. [T][a). The 
vertical arrows mark the position of the boundary at x — 
0. In the range —oo < Xc < 0, note the development of 
exponentially growing tails, associated with the fact that the 
value of V above is very close, but not equal, to an integer 
(here 2). Apart from the tails, both orbitals portrayed here 
are very close to the eigenfunctions of a ID harmonic oscillator 
centered at 3j ^ sec Eq. ©. 



for an unconfined harmonic oscillator; they also exhibit n 
zeros. When i/ n > 0, _D^(+oo) — 0, and for ^ < an 
additional zero (compared to the case oi v = n) develops, 
through which the parabolic cylinder function crosses the 
X-axis and then develops an exponentially growing tail; 
an example for a state with (energy) v — 2-1-9.229 x 10~^^ 
and Xc = 7 [see Fig. [I^a)] is given in Fig. Specifically 
the number of zeros of £'i/(C) (with v > 0) is given by 
the ceiling function^^^ {i^l (this includes the case when 
is a positive integer n); for v < 0, Dy{S^) has no zeros. 

In the case of a semi-infinite graphene sheet with zigzag 
edges (extending for < a: < c»), one can require that 
the additional zero for the xb spinor component coincides 
with the origin of axes (x = 0); this provides the following 
transcendental equation for determining the energy levels 
of the Dirac electrons (remember that v = £^/2): 
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FIG. 3. A double-edge state associated with the third (n = 
2) LL: The xb (top frame) and xa (bottom frame) Dirac- 
spinor components are displayed for Xc = —1 and v = 7.5266 
[see the transcendental Eq. ((8])]. This state lies on the rising 
branch extending out from the flat segment of the third LL 
curve in Fig.[Ua). The vertical dashed lines mark the position 
of the physical boundary at a; = 0. In the range — oo < 
Xc < 0, note the development of exponentially growing tails, 
associated with the fact that the value of i/ above is not equal 
to an integer. In the range — oo < a; < -foo, the number of 
zeros associated with the B component is [v] = 8; for the A 
component it is [r^ — 1] — 7. The physically relevant range 
< X < -l-oo contains only two zeros for both the cases of the 
B and A spinor components. 



The single-particle energies e [which are solutions of 
Eq. ([8])] as a function of Xc are displayed in Fig.Ilfa). One 
sees that Landau levels (with energy e ~ A/2n, n — 0, 1, 
2, . . .) are formed when the centroid Xc of the orbitals is 
far away from the physical edge; the approximate sym- 
bol « signifies that does not take integer values, but 
comes extremely close to them. An illustrative case of 
the corresponding Dirac-spinor orbitals xb and xa are 
portrayed in Fig. [2j In the domain — oo < a; < 0, a tail 
develops due to the fact that the index v is not an in- 
teger. In the physically relevant domain < a: < -l-oo, 
the two components are bulk-like and very similar to the 
familiar wave functions of a ID unconstrained harmonic 
oscillator. Similar orbitals (differing only in the number 
of zeros) apply for all Landau levels with v ^ n > 1. 
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FIG. 4. A bulk-edge state in the LLL (n = 0): The xs (top 
frame) and xa (bottom frame) Dirac-spinor components are 
displayed for Xc = 7 and u = 2.049 x 10~^^ [see the tran- 
scendental Eq. (|8])]. This state lies well inside the quasiflat 
segment of the LLL curve in Fig. [TJa). The vertical arrows 
mark the position of the boundary at j: = 0. In the top frame, 
note the development (in the range — oo < Xc < 0) of an expo- 
nentially growing tail, associated with the fact that the value 
of v above is very close, but not equal, to an integer (here 0). 
Apart from the tail, the orbital portrayed in the top frame is 
very close to the ground-state eigenfunction of a ID harmonic 
oscillator centered at Eq. 0. 



For positive values of Xc near the boundary, and also 
for negative values of Xc, double-edge states are formed 
reminiscent of the single-edge states familiar from the 
theory of the integer quantum Hall effect <^ The corre- 
sponding orbitals for an illustrative case (with = — 1 
and ly = 7.5266) are displayed in Fig. [3] Again, one sees 
the development of a tail in the domain — oo < a: < 0, 
due to the fact that the index v is not an integer. 

The case of large and positive Xc in the LLL (n = 
0) is special and of particular significance regarding 
the strongly correlated Dirac-electron states in finite 
graphene samples under high magnetic field. Indeed in 
this case, the Dirac spinor contains orbitals of both bulk 
and edge character. An illustrative case (with Xc = 7 and 
v = 2.049 X 10~^^ is portayed in Fig.Hl One sees that the 
B component is bulk-like and similar to the ground-state 
of an ID harmonic oscillator in the physically relevant 




FIG. 5. The TB energies at high magnetic field (</> = 0.01) 
for a triangular flake with zigzag edges comprising TV — 6397 
carbon atoms. The n — 0, ±1, ±2, and ±3 Landau levels 
correspond to the flat segments of the curve. The double-edge 
states correspond to the rising-in-energy branches between 
the flat segments. The integer index p counts the TB states 
(negative p values correspond to negative energies) . 



domain < x < +oo. However, in the same domain, the 
A component is clearly edge-like. 

Further understanding of this LLL behavior can be 
achieved through the observation that for Xc >> the 
LLL Xa component can be approximated by 



^LLL,app(^) 



= Ce^'^^~^'=^^y'^erfc(a; - Xc). 



(9) 



Taking into consideration that erfc(— Xc) — 2 for (large) 
Xc » 1, and keeping the lowest order in x/xc in the 
exponent, one can determine the normalization constant 
C. The final simplified expression is: 



(10) 



Eq. (ITUl) has the form of an exponential function de- 
caying inside the graphene sheet. This form agrees very 
well with the full solution of xa in Fig. U [see Eq. ([5]) 
with Xc = 7 and ly = 2.049 x 10~^^]. We note that the 
surface character of this edge-like LLL A-component be- 
comes more pronounced (i.e., it exhibits a narrower 1/xc 
width) the larger the (positive) value of the centroid Xc- 

Solutions for e = 0. In this case, the two equations in 
Eq. ID) decouple yielding the two solutions 



XAiO = 0, 



and 
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(11a) 
(lib) 



(12a) 
(12b) 
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For the nontrivial case {Ca or Cb 7^ 0), neither 
of these two solutions satisfy the boundary conditions 
Xb{-V^Xc) = and xs(+oo) = xa{+oo) = 0. Thus 
there is no dispersionless solution with e = associated 
with the K valley. However, as we will see below, such 
e = solutions exist for the K' valley. 



2. K' valley 




N^6397 
ct)^0.01 




E=0.1699e-07 



The continuous Dirac-Weyl equation coupling the 
and x'si^) components in graphene's K' valley 
is given by 



^Xs - (a;-a:c)XB = (13a) 



^(-Xa) + ix- Xc)i-x'A) = eXs- 



(13b) 



We note that Eq. (|13p has the same form as Eq. ([T]) 
with the substitution xa ^ x'b ^^'^ Xb ^ ~x'a- ^ 
result, for £ ^ 0, one has the following solutions for the 
Dirac spinor in the K' valley 



x'a{0 = ~CMO, 



(14a) 
(14b) 



with the index 1/ — as was the case in the K valley. 

The transcendental equation in the K' valley for the 
indices v (or energies e) as a function of Xc is written as 



(15) 



The solutions of Eq. (|15p as a function of plotted 
in Fig. lUb). Note that the index v > 1 m this case 
[-D^_i(^) has no zeros for ly' < I]. This contrasts with the 
case of the K valley shown in Fig. [Tfa), where z/ > 0. 

For e = 0, the two equations in Eq. ((T3)) decouple, and 
there is a physically valid solution 



XAiO 

x'siO 



0. 



(16a) 
(16b) 



Assuming a relation i/ = this dispersionless band of 
edge states can be associated with an index v = 0, and it 
is denoted by a thick dashed line in Fig.[IJb). This band 
maintains also for zero-magnetic field, since 



lim e ^ oc e 



(17) 



which represents an edge state for qy < 0; the tilded 
X denotes the x position in the original dimensions of 
length (before the introduction of the reduced variable 
X = x/Ib] see Ref. 0)- 




FIG. 6. Examples of TB electron densities at high magnetic 
field ((/> = 0.01) for the A (left) and B (right) sublattices 
associated with: (a) + (b) mixed bulk-edge dispersive states in 
the n = Q Landau level (whose TB energies reside on the flat 
step at _E ~ in Fig. [5]). (c) double-edge states situated on the 
rising-in-energy branch between the n = and n — 1 Landau 
levels. Energies in units of the hopping coupling parameter t. 



B. Tight-binding approach for finite graphene 
flakes 



In the tight-binding (TB) calculations, we use the 
hamiltonian 



H' 



TB 



(18) 



<hj> 



with <> indicating summation over the nearest- 
neighbor— sites The hopping matrix element 

kj =texp(^j J^' ds- A(r)^ , (19) 

where t — 2.7 eV,ri and rj are the positions of the carbon 
atoms i and j, respectively, and A is the vector potential 
associated with the applied perpendicular magnetic field 
B. 



The calculations were carried out for two shapes that 
support a zigzag edge on all sides of the graphene flake, 
that is, equilateral triangles and regular hexagons, as well 
as for a square shape which exhibits both zigzag and 
armchair edges. The number of carbon atoms consid- 
ered is = 6397 for the triangular flakes, N ~ 6144 
for the hexagonal ones, and N = 2074 for the square 
one. The diagonalization of the TB hamiltonian [Eq. 
(1181) ] is implemented with the use of the sparse-matrix 
solver ARPACK,^ 



III. RESULTS OF TIGHT-BINDING 
CALCULATIONS AND THEIR 
INTERPRETATION 

A. Trigonal graphene flakes 

In Fig. \E\ we display the TB energies for the trian- 
gular flake with (j) = 0.01, (j) = eSB/{hc) being the di- 
mensionless magnetic flux through an hexagonal unit of 
the two-dimensional graphene lattice; S is the area en- 
closed by the hexagon and B is the magnetic field. The 
value (j) = 0.01 corresponds to a magnetic field sufficiently 
high so that Landau levels have been formed, but at 
the same time low enough so that Hofstadter-butterfly^ 
effects (due to the periodicity of the lattice) have not 
developed)^ The TB-energy curve in Fig. [S] exhibits a 
well defined trend, i.e., several almost- flat horizontal seg- 
ments are connected via fast- varying and rising-in-energy 
branches. The flat segments correspond to the Landau 
levels with energy E/t cx sign(n) i/jnl , n = 0, ±1, ±2, 

A close inspection of the properties and behavior of 
the states associated with the TB energies in the n — 
level reveals that this level contains two different bands: 
(i) a nondispersive one with energies close to the avail- 
able machine precision {E/t ^ lO^^'*) [that correspond 
to the e = solutions of the continuous model; see Sec. 
Ill Aj and (ii) a quasidegenerate dispersive band with en- 
ergies that are still very small (starting at E/t ^ 10^^) 
but increase gradually and then merge with the rising- 
in-energy branch; this band corresponds [see the tran- 
scendental equation ([5])] to the e ^ solutions of the 
continuous model in the flat region of the lowest curve in 
Fig. Ufa). The nondispersive zero states coincide with the 
midgap surface states under fleld-free conditions and they 
have been studied extensively)^"— these states are not 
influenced by the magnetic field and will not be discussed 
any further here. The close-to-zero dispersive states are 
formed as a consequence of the external magnetic field, 
and they will be the primary focus of this paper. 

In Fig. [HI we display the TB electron densities (specif- 
ically the square root of the densities"^^) for states asso- 
ciated with the LLL in the case of a triangular graphene 
flake under the same magnetic fleld 4> = 0.01. The elec- 
tron densities of the A (left, online red) and B (right, 
online blue) sublattices are plotted separately. Figs. 
[SJa) and EUb) correspond to TB states with energies 
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FIG. 7. Examples of TB electron densities at high magnetic 
field {cj> = 0.01) for the A (left) and B (right) sublattices 
associated with: (a) + (b) double- bulk states in the n = 1 
Landau level, (c) double-edge states situated on the rising- 
in-energy branch between the n — 1 and n — 2 Landau levels. 
Energies in units of the hopping coupling parameter t. 



{E/t = 0.1699 X lO"'^ and E/t = 0.1475 x 10"^^ re- 
spectively) that lie well inside the flat segment of the 
LLL {n = 0) energy curve in Fig. [51 It is apparent that 
these states are of a mixed bulk-edge character, with 
the B-sublattice component being bulk-like and the A- 
sublattice component being edge-like. They are anal- 
ogous to the mixed bulk-edge LLL states described in 
Sec. Ill Al within the continuous relativistic Dirac-Weyl- 
equation framework. In particular, the spatial profile of 
the A-sublattice component of the TB densities in Figs. 
[Sfa) and [IKb) agrees well with the surface-state Dirac- 
spinor component in Eq. ([TU)) : see also Fig. |31 bottom 
frame. Moreover, the profiles of the TB densities for the 
B sublattices in these figures exhibit the qualitative be- 
havior of a D^[\/2{x — Xc)] function with described 
in Sec. Ill A[ see also Fig. SI top frame. Note that a lower 
TB energy [case of Fig. [D^a)] corresponds to a continu- 
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FIG. 8. Examples of TB electron densities at high magnetic 
field {(f) = 0.01) for the A (left) and B (right) sublattices of 
a triangular flake with A'^ — 6397 carbon atoms and zigzag 
edges, (a) double-bulk states in the n — 2 Landau level, 
(b) double-edge states situated on the rising-in-energy branch 
between the n = 2 and n = 3 Landau levels. Energies in units 
of the hopping coupling parameter t. 



ous D^[^/2{x — Xc)] state with a larger centroid Xc > 0, 
farther away from the physical edge. For states with 
higher TB energies, lying on the rising-in-energy branch, 
the B-sublattice component moves towards the physical 
edge and transforms into an edge state, as illustrated 
by the double-edge TB state in Fig. IHlJc) (with energy 
E/t = 0.3033). 

In Figs. Elja) and[7{b), we portray TB densities for 
two states associated with the second Landau level (with 
index n = 1 at an energy E/t 0.32685). Fig. [Tja) 
corresponds to a Dirac spinor in the K valley having an 
A component consisting of a -DolO state (no nodes in- 
side the graphene flake) and a B component consisting of 
a Di{£_) state (a single node inside the graphene flake). 
Fig. Efb) portrays a similar TB state in the K' valley, 
since the A and B sublattices correspond to the contin- 
uous functions -Di(C) and -Do(C)i respectively, which is 
the opposite from the i^T- valley case in Fig.[71[a); see Sec. 
Ill A 21 Fig.[71Jc) portrays a double-edge state with energy 
E/t — 0.384558 (lying in Fig. [5] on the rising-in-energy 
branch between the n — 1 and n — 2 Landau levels). 
It is apparent that this double-edge TB state is associ- 
ated with the K' valley and has been evolved out of the 
double-bulk state in Fig. El^b); note the preservation of 
the single node (no node) topology in the A-sublattice 
(B-sublattice) electron-density component. 

In Fig. [SJa) , we portray TB densities for a double-bulk 
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FIG. 9. The TB energies at high magnetic held (0 = 0.01) 
for an hexagonal graphene flake with zigzag edges comprising 
N = 6144 carbon atoms. The n = 0, ±1, and ±2 Landau lev- 
els correspond to the flat segments of the curve. The double- 
edge states correspond to the rising-in-energy branches be- 
tween the flat segments. The integer index p counts the TB 
states (negative p values correspond to negative energies) . 

state associated with the third Landau level (with index 
n = 2 at an energy E/t « 0.458268; see Fig. E]). The TB 
densities in Fig. [DJa) correspond to a continuous Dirac 
spinor in the K' valley having an A component consist- 
ing of a Z?2 (0 state (two nodes inside the graphene flake) 
and a B component consisting of a Di{S^) state (a single 
node inside the graphene flake); see Sec. Ill A 21 and Fig. 
Fig. [2] [but with B replaced by A' and A replaced by 
B' (online: blue -n- red)]. Fig. Mjo) portrays a double- 
edge state with energy E/t = 0.501967 (lying in Fig. [5] 
on the rising-in-energy branch between the n = 2 and 
n = 3 Landau levels). It is apparent that this double- 
edge TB state is associated also with the K' valley and 
has been evolved out of the double-bulk state in Fig.jSja); 
note the preservation of the two-node (single- node) topol- 
ogy in the A-sublattice (B-sublattice) electron-density 
component. Note the similarities with the double-edge 
continuous Dirac spinor within the physical semi-infinite 
graphene plane portrayed in Fig. |3l 



B. Hexagonal graphene flakes 

In Sec. nil Al we studied the nature of the TB states 
in the case of a trigonal flake with zigzag edge termi- 
nations. A characteristic property of triangular flakes is 
that the same sublattice participates in the edge termina- 
tions of two adjacent polygonal sides (forming an angle 
of 60°). In this section, we study hexagonal graphene 
flakes^^ with zigzag edge terminations, which is a more 
complicated case. This is due to the fact that the sublat- 
tices A and B alternate in providing the edge termination 
of adjacent polygonal sides (having an angle of 120°). 

In Fig. [9l we display the TB energies for an hexago- 
nal flake with (p — 0.01 (corresponding to a high mag- 
netic fleld where Landau levels have been formed). As 
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was the case with the trigonal TB energies in Fig. [5l 
the TB-energy curve in Fig. [3] exhibits also a well de- 
fined trend, i.e., several almost-flat horizontal segments 
which are connected via fast-varying (rising-in-energy) 
branches. The flat segments correspond to the Landau 
levels with energies E oc sign{n)y^\n\, n — 0, ±1, ±2, . . .. 

Furthermore, as was also the case with the trigonal 
flakes, a close inspection of the properties and behavior 
of the states associated with the TB energies in the n ^ Q 
level in Fig. [S] reveals that this level contains two different 
bands: (i) a nondispersive one with energies close to the 
available machine precision {E/t ~ 10"^'') [that corre- 
spond to the e = solutions of the continuous model; see 
Sec. Ill A] and (ii) a quasidegenerate dispersive band with 
energies that are still very small (starting at E/t ^ 10~^) 
but increase gradually and then merge with the rising en- 
ergy branch; this band corresponds [see the transcenden- 
tal equation ([5])] to the £ ^ solutions of the continuous 
model in the flat region of the lowest curve (i.e., in the 
LLL) in Fig. [T](a). Here, we study these dispersive states 
in the LLL that are formed due to the presence of the 
magnetic fleld. 

In Fig.[TUl we display the TB electron densities (specif- 
ically the square root of the densities—) for states associ- 
ated with the LLL in the case of an hexagonal graphene 
flake under the same magnetic fleld 4> = 0.01. The elec- 
tron densities of the A (left, online red) and B (right, 
online blue) sublattices are plotted separately. Fig.fTOla) 
corresponds to a TB state with energy {E/t = 0.1459 x 
10"'^) that lies well inside the flat segment of the LLL 
(n = 0) energy curve in Fig. [S] This state is of a mixed 
bulk-edge character, exhibiting, however, a more complex 
proflle compared to the corresponding mixed bulk-edge 
states for the trigonal flake in Figs. IHl^a) andlHKb). This 
is due to the alternation of the A and B sublattices along 
the polygonal sides forming the edge of the hexagonal 
flake. In particular, focusing on a side with a B-sublattice 
termination (e.g., the one at the upper-right corner), one 
sees that the A-sublattice density (left, online red) ex- 
hibits an edge-state behavior, while the B-sublattice den- 
sity (right, online blue) exhibits a bulk-state proflle. Fo- 
cusing on a side with an A-sublattice termination (e.g., 
the one at the lower- left corner), one sees the opposite, 
i.e., the A-sublattice density (left, online red) exhibits a 
bulk-state behavior, while the B-sublattice density (right, 
online blue) exhibits an edge-state proflle. In the continu- 
ous relativistic Dirac-Weyl model (Sec. Ill A[) . the former 
case is associated with spinor components in the K val- 
ley (permitting the vanishing of xb on the edge; see Sec. 
Ill A ip . while the latter is associated with spinor compo- 
nents in the K' valley (permitting the vanishing of x'a 
on the edge; see Sec. Ill A2[) . 

Fig. [TUf b) displays the TB densities for states with a 
higher TB energy, E/t — 0.2576, lying on the rising-in- 
energy branch between the n = and n — 1 Landau 
levels in Fig. [HI This state represents a double-edged one 
and it can be interpreted as having continuously evolved 
form the state portrayed in Fig. [TUf al. with the bulk 
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E=0.1459e-03 




FIG. 10. Examples of TB electron densities at high magnetic 
field {(j) = 0.01) for the A (left) and B (right) sublattices of 
an hexagonal flake with A'^ — 6144 carbon atoms and zigzag 
edges, (a) a mixed bulk-edge state in the n — Landau level, 
(b) a double-edge state with energy situated on the rising-in- 
energy branch between the n = and n = 1 Landau levels. 
The labels A and B indicate the (alternating) sublattice-type 
edge termination along the sides forming the physical bound- 
ary of the hexagon. Energies in units of the hopping coupling 
parameter t. 



components having been pushed against the edges. 

Fig. [11] portrays TB densities for two double-bulk 
states associated with the second Landau level of the 
hexagonal flake (with index n = 1 at an energy E/t ^ 
0.3270; see Fig. The sublattice densities in Fig. 

[TTT a) correspond to a TB state with the lowest energy 
{E/t K, 0.3268) in this Landau level, and they are concen- 
trated around the center of the flake, far away from the 
edges. In addition, neither one of them (A-sublattice den- 
sity or B-sublattice density) exhibits any nodes, which 
seems paradoxical at a first glance for a state belonging 
to the first Landau level. The explanation can be found 
in that the edge of the flake has very little influence in 
this case, which thus corresponds to a nodeless I = 
Dirac-Weyl state of the n = 1 LL in a circular graphene 
dot with zigzag termination; see Eqs. (A. 2) and (A. 4) in 
Ref . 113. 

Fig. HHb) portrays TB densities for another double- 
bulk state with energy E/t = 0.3271 belonging to the 
second (n = 1) Landau level of the hexagonal flake; see 
Fig. [HI The density proflles in Fig. fTTTb) exhibit a zero- 










E=0.3271 




E=0.3323 




FIG. 11. Examples of TB electron densities at high magnetic 
field {(/} = 0.01) for the A (left) and B (right) sublattices of 
an hexagonal flake with N = 6144 carbon atoms and zigzag 
edges. Both rows [(a) and (b)] portray double-bulk spinors 
in the n = 1 Landau level. The labels A and B indicate the 
(alternating) sublattice-type edge termination along the sides 
forming the physical boundary of the hexagon. Energies in 
units of the hopping coupling parameter t. 



FIG. 12. Examples of TB electron densities at high magnetic 
field {(f) = 0.01) for the A (left) and B (right) sublattices of 
an hexagonal flake with A'^ = 6144 carbon atoms and zigzag 
edges. Both rows [(a) and (b)] portray double-edge spinors 
with energies in the rising-in-energy branch between the n = 1 
and n = 2 Landau level. The labels A and B indicate the 
(alternating) sublattice-type edge termination along the sides 
forming the physical boundary of the hexagon. Energies in 
units of the hopping coupling parameter t. 



node and a one- node structure in analogy with the -Do(0 
and -Di(^) functions describing the n = 1 Landau level in 
the continuous Dirac-Weyl model (see Sec. Ill Al) . How- 
ever, in contrast to the single-edge semiinfinite graphene 
plane in Sec. Ill Al in Fig. fTTTb) the zero- node/one- node 
topology is present in both sublattices (in both the left 
and right panels) due to the alternation of the edge ter- 
mination along the hexagon's sides between the A and B 
sublattices. 

Fig. [T2] portrays TB densities for two double-edge 
states whose energies lie on the rising-in-energy branch 
between the n = 1 and n = 2 Landau levels; see Fig. [HI 
The sublattice densities in Fig.[T2la) correspond to a TB 
state with energy E/t = 0.3756. It is apparent that they 
can be viewed as having evolved out of the densities in 
Fig. mT b). with the centroids of the bulk densities hav- 
ing been pushed against the hexagonal edges. Remark- 
ably, at the same time, the zero- node/one- node alternat- 
ing nodal topology is preserved in the double-edge state 
in Fig. IT^ a) in complete ananlogy with the double-bulk 
state in Fig.lTTTb). 

The sublattice densities in Fig. IT^ b) correspond to a 
TB state with energy E/t = 0.3323; they also represent 



a double-edge state like the one in Fig. IT^a). However, 
a new characteristic [compared to the double-edge state 
in Fig. IT^ a)] is the appearance of an additional nodal 
pattern resulting from wave-function quantization along 
(and parallel to) the hexagon's sides. This additional pat- 
tern is superimposed upon the perpendicular-to-the-edge 
zero- node/one- node pattern (the latter being present al- 
ready in the continuous Dirac-Weyl model of a semiinfi- 
nite graphene plane studied in Sec. Ill A[) . We note that 
such combined parallel-to and perpendicular-to-the-edge 
nodal patterns (not shown) have been found by us also 
in several instances of TB states in the case of a trigonal 
graphene flake in Sec. IIII Al 



C. Square graphene flakes 

In this section, we study the more complicated case 
of a square graphene flake which necessarily has mixed 
armchair and zigzag edge terminations (each type of ter- 
mination developing on opposite sides of the square). In 
Fig.lHl we display the TB energies for such a square fiake 
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FIG. 13. The TB energies at high magnetic field {(p = 0.01) 
for a square graphene flake with both zigzag and armchair 
edges comprising A'^ — 2074 carbon atoms. The lowest (n — 0) 
Landau level corresponds to the fiat segment of the curve at 
i? ~ 0. The integer index p counts the TB states (negative p 
values correspond to negative energies). 



with = 0.01 as in the previous studied cases of trigonal 
(Sec. IIII Ap and hexagonal flakes (Sec. IIII B|) . Compared 
to the energy curves in the previous cases [trigonal (see 
Fig. [SJ and hexagonal (see Fig. [3])], the TB energies in 
Fig. [T3l exhibit higher Landau levels (horizontal segments 
in Fig. [T3| with n > 1 that are not as well formed; this 
may be due to the smaller number of carbon atoms in the 
square flake {N = 2074). The LLL (n = 0) level, how- 
ever, is well formed, and this is sufficient for our purposes 
here, namely to investigate whether the mixed LLL bulk- 
edge states maintain in the presence of edge segments 
with armchair termination. 

Indeed the TB sublattice densities for the LLL state 
(with energy E/t = 0.2769 x 10"^) portrayed in Fig. [H 
show that the mixed LLL bulk-edge behavior maintains 
also in the case of a square flake. Naturally, due to the 
coupling between the K and K' valleys induced by the 
presence of the armchair terminations, each sublattice 
[A (online red) and B (online blue)] exhibits now both 
edge and bulk density contributions, which however cor- 
respond to different valleys as explicitly marked in Fig. 

M 

Mixed bulk-edge states of a square graphene flake in a 
perpendicular magnetic field were also reported in a re- 
cent study. In this study the appearance of such mixed 
states with significant weight at the zigzag edges were 
attributed to the coupling between the K and K' valleys 
due to the armchair edges (see in particular Sec. IV in 
Ref . [26h . This interpretation differs from the conclusion 
presented in our paper where the occurrence of mixed 
bulk-edge LLL states is shown to originate solely from 
the zigzag edge termination. 



N=2074 (t)=0.01 




E=0.2769e-02 

FIG. 14. Example of TB electron densities at high magnetic 
field (0 = 0.01) for the A (online red; top frame) and B (on- 
line blue; bottom frame) sublattices of a square flake with 
A'' = 2074 carbon atoms and both zigzag and armchair edges. 
A state representing a mixed bulk-edge state in the LLL is 
displayed. The density contributions on each sublattice asso- 
ciated with a given valley are also marked. Energies in units 
of the hopping coupling parameter t. 



IV. SUMMARY AND DISCUSSION 

The properties of single-electron states in graphene 
flakes with zigzag edge termination under high magnetic 
fields (in the regime of Landau-level formation) were in- 
vestigated using tight-binding calculations. A systematic 
interpretation of their character (bulk-like versus edge- 
like) was achieved via a comparison of the tight-binding 
electron densities with analytic expressions (based on 
parabolic cylinder functions) for the relativistic Dirac- 
Weyl spinors in the case of a semi-infinite graphene plane. 
A variery of graphene flakes was considered, namely, trig- 
onal, hexagonal, and square ones. 

The higher Landau levels were found to comprise ex- 
clusively electrons of bulk-type character (for both sub- 
lattices). Furthermore, electrons with energies on the 
rising-in-energy branches (connecting the Landau levels) 
are described by edge-type states reminiscent of those fa- 
miliar from the theory of the integer quantum Hall effect 
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for nonrelativistic electrons. 

In contrast, in all cases studied the lowest, {n ~ 0) 
Landau level contained relativistic Dirac electrons of a 
mixed bulk-edge character without an analog in the non- 
relativistic case. Most importantly, it was shown that 
such mixed bulk-edge states maintain also in the case of 
a square flake with combined zigzag and armchair edge 
terminations. 

The presence of mixed bulk-edge LLL states in 
graphene samples with realistic shapes points at sig- 
nificant implications concerning the many-body corre- 
lated FQHE excitations. We recall that Ref. 17 studied 
the many-body correlated FQHE excitations in the LLL 
in the simplified case of a circular graphene flake with 
zigzag edge termination, and it found that the two-body 
Coulomb-interaction matrix elements are given as a sum 
of four terms 

^((^i^2|fo3^4) + (eie2|e3e4) + (616216364) + (616216364)), 

(20) 

where |6) and |e) denote the bulk and edge components of 
the mixed LLL state; note that, due to the equal weights 
(50% — 50%) of the bulk-like and edge-like components, a 
prefactor of 1 /4 appears in front of each term in Eq. ([20| . 
As a consequence of Eq. (PH)) and of the 1/4 prefactor, a 
sizable attenuation of the many-body correlated FQHE 



excitations in the LLL (associated with the (6162I6364) /4 
term reflecting the depletion of the bulk component) was 
found in the simplified case of a circular graphene fiake. 
Furthermore, it was shown^^ that the insulating behavior 
at the Dirac neutrality point under high B (experimen- 
tally observediSiiS in graphene samples along with the 
1/3 FQHE) is associated with the Coulombic repulsion 
due to the accumulation of charge at the edges [related 
to the (6162 1 6364) /4 and the remaining two cross terms 
in the Coulomb-interaction matrix elements given in Eq. 

mi 

The current study shows that the appearance of mixed 
bulk-edge LLL states is a property of the presence of seg- 
ments in the graphene-sample boundary having a zigzag 
edge termination, independent of the precise shape of the 
graphene sample. This finding suggests that the results 
of Ref. concerning the many-body correlated FQHE 
excitations in the LLL can be generalized to graphene 
samples with more realistic shapes. 
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